Renormalization group method for weakly-coupled quantum chains: application to the 

spin one-half Heisenberg model 
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The Kato-Bloch perturbation formalism is used to present a density-matrix renormalization-group 
(DMRG) method for strongly anisotropic two-dimensional systems. This method is used to study 
Heisenberg chains weakly coupled by the transverse couplings J± and Jd ( along the diagonals) . An 
extensive comparison of the renormalization group and quantum Monte Carlo results for parameters 
where the simulations by the latter method are possible shows a very good agreement between the 
two methods. It is found, by analyzing ground state energies and spin-spin correlation functions, 
that there is a transition between two ordered magnetic states. When Jd/J± ^ 0.5, the ground state 
displays a Neel order. When Jd/J± ^ 0.5, a collinear magnetic ground state in which interchain spin 
correlations are ferromagnetic becomes stable. In the vicinity of the transition point, Jd/J± ~ 0.5, 
the ground state is disordered. But, the nature of this disordered ground state is unclear. While the 
numerical data seem to show that the chains are disconnected, the possibility of a genuine disordered 
two-dimensional state, hidden by finite size effects, cannot be excluded. 



I. INTRODUCTION 

In a recent publication!, it was shown that the density- 
matrix renormalization group method (DMRG)^ can be 
applied to an array of weakly coupled quantum chains. 
As an illustration of the method, weakly coupled Heisen- 
berg spin chains were studied and some partial results on 
the ground state energies were shown to be in good agree- 
ment with previous quantum Monte Carlo (QMC) stud- 
ies. But the essential question concerning the stability of 
the disordered one-dimensional (ID) ground state against 
small transverse perturbations was not addressed. 

The motivation behind such a study is in the search of 
a disordered ground state for a spin one-half Heisenberg 
model in dimension higher than one. A spin liquid state 
without spin rotational or translational symmetry break- 
ing has been conjectured to be relevant for the physics 
of high-Tc cuprate superconductors^. A possible candi- 
date is the resonance valence bond (RVB) state 5 . Earlier 
attempts^^* 8 ^ to find the RVB ground state by various 
techniques (l/S expansions, exact diagonalization, quan- 
tum Monte Carlo) have given some indication about its 
possible realization. But their conclusions are still dis- 
puted. It has even been arguedi - that the spin-Peierls 
mechanism, not RVB, may be the most natural way to 
lead to a disordered state. 

More recently, the interest has shifted to search for 
a RVB state on quasi ID systems. A pure spin one- 
half Heisenberg chain has a disordered ground state with 
neutral spin one-half excitations (spinons) and does not 
break spin rotational or translational symmetry. It is 
thus tempting to try to find a higher dimensional gener- 
alization of this state by the application of small pertur- 
bations. Contrary to an earlier claim of the realization 
of a spin liquid 11 , subsequent studiesi2iiii4 indicate that 
the introduction of the rung transverse coupling J± be- 
tween the chains (see Fig.^) seems to lead to a Neel state 
for any non-zero J± . A possible way to avoid the Neel or- 
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FIG. 1: Sketch of the ground state of weakly coupled Heisen- 
berg chains as function of J± (dashed lines along the rungs) 
and Jd (dotted lines along the diagonals): Neel state when 
Jd/J± 0-5 (left), collinear state when Jd/J± ^ 0.5 (right) 



der is to introduce, in addition to J±, a small frustration 
Jd along the diagonals. In a recent worki 5 . it was claimed 
that a spin liquid state is realized when J± = 2Jd- 

A more direct motivation in studying a model of 
weakly coupled Heisenberg chains stems to its relevance 
to the understanding of interchain effects in quasi-one- 
dimensional materials 16,17 ' 18 . A recent neutron scat- 
tering cxpcrimenljiS on the frustrated antiferromagnet 
(AFM) Cs^CuCl^ found that the dynamical correlation 
show a highly dispersive continuum of a excitations with 
fractional quantum numbers, a signature of a spin liquid 
state. 

In this paper, a general formalism of the DMRG algo- 
rithm for weakly coupled chains of Ref. 1 is presented. 
This method is a particular case of a recent matrix 
version^ , of the general perturbation expansion which 
was proposed decades ago by Kato and BlochSii 2 ^* 2 ^. 
The Kato-Bloch expansion was initially introduced to 
find the correction on a single state. This expansion is 
straightforwardly generalized to account for many low 
lying states. The method is in spirit close to an ear- 
lier perturbative renormalization group by Hirsch and 
Mazenke 24 . A more detailed study of weakly coupled 
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Heisenberg chains is presented. An extensive compari- 
son with quantum Monte Carlo results for unfrustrated 
transverse couplings is made. It shows a good agree- 
ment between the two methods when the perturbation is 
small and the lattice not too large. Then the question of 
the stability of the non-magnetic state in the presence of 
frustration is addressed. It is shown that the 2D DMRG 
algorithm can provide a convincing answer to this ques- 
tion, at least for the parameters that were investigated. 
It is found, by analyzing ground state energies and spin- 
spin correlation functions, that the perturbation is rel- 
evant, leading to magnetic ground states (see Fig. 
When Jd/J± ^ 0.5, the ground state displays a Neel 
order. When Jd/J± ^ 0.5 a collinear magnetic ground 
state in which interchain spin correlations are ferromag- 
netic becomes stable. In the vicinity of the transition 
point, Jd/J_L ~ 0.5, the system seems to behave as an 
assembly of independent chains. This is reminescent of 
the so-called sliding Luttinger liquid 26 recently found in a 
model of crossed spin one-half Heisenberg chains^. But 
it is impossible to exclude a genuine 2D spin liquid state 
(i.e., with a spin gap) masked by finite size effects. 



II. FORMAL DEVELOPMENT 

The DMRG method described below can work for spin, 
fcrmionic as well as bosonic systems, and so it is conve- 
nient to use a general formulation of the algorithm that 
can then be adapted to each of these cases. The model 
Hamiltonians under consideration can be written as fol- 
lows: 



and the set of approximate eigenvalues of Hn is given by 
the sum 



E\\[ n ] = e ni + e n2 + ... + e nL , (4) 

where [n] = {ri\, 712, til) and n\ labels to an eigenset 
on the chain I. 

Let P|| be the projector on the states 

P \\ =X) I %»]><* UN I ( 5 ) 

[n] 

and Q|| = 1 — P\\ . 

Let (Er n ], |3>[ n ])) be the exact eigenset of H. This 
eigenset will tend to (£J||r n i, |$iir n i)) in the limit g — > 0. 
Let P be the projector onto the spaces P may be 

written as follows 

N 

Since the perturbation g is small, it is assumed that the 
subspaces generated by the |$ii[„])'s and by the |$r n i)'s 
are not orthogonal. An approximate expression of HP 
in the basis spanned by the eigenstates of Hn will now 
be derived by using a generalization of a method first in- 
troduced by Kato2i and later modified by Bloch 2 ^. The 
advantage of the Bloch's version is that it leads to a sim- 
pler expansion. 

Following Bloch, let U be the operator 



H = H ll +gH ± , (1) 

where H\\ is a the sum over one-dimensional (ID) Hamil- 
tonians (longitudinal direction), 

L 

#11 =E^ ( 2 ) 

i=i 

and H± is the interaction between these ID systems 
(transverse direction). The coupling constant g is such 
that g < 1. 

Since g -C 1, it is natural to study the problem using 
perturbation theory. The Kato-Bloch formalism is con- 
venient to set up a perturbation expansion around a nu- 
merical solution of Hn provided by the DMRG. For a sin- 
gle chain I whose Hamiltonian is Hi , a set of eigenstates 
\4>m) an d eigenvalues e ni can be obtained by the usual 
ID DMRG. The zeroth order set of eigenstates I^iim) of 
the full longitudinal Hamiltonian is simply given by the 
tensor product of the \<j> ni ), 

I*iin) = |0 ni )|0„ 2 )...|0 ni ), (3) 



N 

which projects the l^iirnl) onto |$r n i), and U satifies, 
UP\\ — hi. The problem of finding an expansion of HP 
projected onto Pii is equivalent to finding an expansion 
iovP\\HU. 

One starts by deriving an equation satisfied by U. The 
Schrodinger equation 

JJ|$ [n] > =£[„]!$[„]> (8) 
is transformed as follows, 

{H - H^ [n] ) = £ [n] \^ [n] ) (9) 

where H\\ is identical to Hn in the subspace spanned by 
thc l$||[n])'s, and, 

£ [n] =E ln] -(^ [n] \H^ [n] ). (10) 
When a single state |$||[o]) is kept, Hn is given by 
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This leads to the expansion for U 
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The method reduces to the usual stationary perturbation 
expansion. It is known that such an expansion does not 
often converge. The main source of divergence is the near 
degeneracy of the eigenvalues. Now if many states up to 
a cut-off n c are kept, a possible generalization of Hu to 



many states |3>|i[o] > ■ ■ ■ \<& 
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Thus if n c is suitably chosen, the series will converge^. 
The purpose of this choice is to shield the eigenvalue E^ 
from the rest of the spectrum by treating the n c — 1 states 
just above the ground state exactly and the remaining 
spectrum perturbatively. 

By applying Pit and then U to the Equation([5J| above, 
one finds, 



gUH ± \^. 



= £r».il$r. 



(11) 



The subtraction of Equation( 111(1 from Equation( |SJ) 
leads to 



{H-Hu -gUH ± )\^,. 



0. 



(12) 



By applying (^nfri]! on the right of eauationf I12[) and 
performing the summation over [n] , one finally obtains 
the equation satisfied by U, 



(H - Hu - gUH ± )U = 0. 



(13) 



Eauationf I13fl is further transformed by using the fact 
that P\\U = P|| and U = P\\U + Q\\U. One obtains: 



U = P\\+ gQ\\ {HJA - UHJA) 



where Qu is given by 



Oil =Q||(-ff||-JI||)- 1 . 



(14) 



(15) 



«<°> = P 



n-l 



(16) 



U {n) = gQ\\ [H±U {n - x) - U^HxU^-p-^] (17) 



0=1 



From this expansion, one finds the approximate Hamil- 
tonian H = P\\HU is 



H = ^S|| N |$|| N )<$|| W | + gP ll H x P ll + 

[n] 

g 2 P l] H ± Q ll H ± P ]l +... (18) 

This perturbation expansion is a matrix generalization 
to many states of Bloch's expansionist which was estab- 
lished for a single state. Even though the ground state 
and a few low lying states will ultimately be computed, it 
is important to keep many low lying states in the pertur- 
bation expansion. This is because the convergence will 
mainly depend on two quantities. The first one is obvi- 
ously g. The second one is the projector Q\\ . If Hu = Hu , 
then Q|| = 0. In that case only the first order term in 
equation f USf) is not equal to zero. The rewriting of the 
original problem to equation ( I18|) is a simple change of 
basis. So in the limit where n c = dimA, where A is 
the Hilbert space in which all the operators are defined, 
the method is exact. But since only a small number of 
eigenstates of Hu can be used even if the full spectrum is 
known, Qu =/= 0. The magnitude of Qu in the expansion 
decreases by increasing the cut-off n c . It is to be noted 
this matrix expansion is close to the method of Hirsch 
and Mazenko^*, who also used a block expansion near 
the solution of an unperturbed Hamiltonian. The prob- 
lem with their study was, however, that their technique 
was applied to a model with no small parameter. 

When the DMRG is used as a method of solution for 



Hu, we can not know Q\\ exactly. This is because the 



DMRG does not keep any information about the trun- 
cated states. But it is possible to define a perturbative 
expansion in a reduced space spanned by the states kept. 
The above perturbative expansion will thus be adapted 
in this study as follows. During the ID DMRG part of 
the method, N s — ms\ x msi states will be obtained 
for the reduced superblock (i.e., the superblock reduced 
to the two external large blocks; it is supposed that 
open boundary conditions (OBC) are used). Typically 
msi = 16 — 192 during this investigation. The complete 
spectrum of this reduced superblock can be obtained as 
in the thermodynamic algorithm 27 . This spectrum will 
serve as Hu . Only a small fraction msi = 16 — 96 of 
these states can be kept for the generation of the 2D 
lattice. The ms 2 states will define Pu, and Qu is con- 



structed using the remaining Ns — ms2 states. Hence 
the perturbation expansion in Eq. ( I18|) will be made by 
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assuming that H\\ is the low energy Hamiltonian of size 
msi x msi obtained from the DMRG rather than the 
exact ID solution of H\\ . 

The Hamiltonian H is one-dimensional and it will be 
studied by the DMRG method. The only difference with 
a normal ID situation is that the local operators are now 
ms2 x ms2 matrices which makes computations heav- 
ier, ft should be noted that the accuracy of the method 
is related to the diagonalized unperturbed Hamiltonian 
obtained from the DMRG. This Hamiltonian, although 
it leads to a very accurate ground state energy, is less 
accurate for high lying states and correlation functions. 
So the potential errors of the method will come from 
the DMRG as well as the truncated perturbative se- 
ries. A better approach is to use the exact diagonal- 
ization method to diagonalize the unperturbed Hamilto- 
nian. However, in that case one will be restricted to small 
chains. 



III. APPLICATION TO THE HEISENBERG 
MODEL 

The above formalism will now be applied to the 
anisotropic Heisenberg model on a 2D square lattice. The 
Hamiltonian reads: 

Hspins — S-ySj+l,; + J± SijS^j + i + 

i,l i,l 

Jd / ](Sj,iSt+i,;+i + Si+i^Sj.i+i) (19) 
i,l 

where the Sij are the usual spin one-half operators. 

The question of the condition of the onset of long-range 
order as a function of Jj_ has been addressed in many 
studies. Spin-wave analysi o 11 ! 28 predicted that there is a 
finite critical J± c ~ 0.03, above which long-range order 
is established. Renormalization group analysis^ supple- 
mented by series expansion computations found that if 
J± c is finite, it cannot exceeds 0.02. A finite critical 
value is at variance with a random phase approximation 

I 



The second order term ( I23|) generates a long-range 
coupling between the chains, which makes it difficult to 
treat. One can see that the condition for the matrix ele- 



(RPA) 13 which predits J ±c = 0. The QMC method com- 
bined with a multichain mean-field approach— has con- 
cluded that when J± — 0, the ground state is an anti- 
ferromagnet down to J± = 0.02. From these studies, it 
is likely that the AFM ground state is stable as soon as 
J± 7^ 0. This does not, however, preclude a spin liquid 
ground state in the case when Jd is added between the 
chains. When this exchange term is added, the QMC 
method faces the infamous sign problem. The two-step 
DMRG method presented here can help to find, if it ex- 
ists, the spin liquid ground state. 

The adaptation of the formalism discussed in sec- 
tion^ |Hj to the model of Equation( I19|l is without any 
difficulty. The first step is the solution of the ID Hamil- 
tonian: 

H l = J2Si,iSi+i,i (20) 

i 

by the usual DMRG method. This yields the chain eigen- 
values e„, and eigenstates | </>«,). From equation f lTBj) . the 
projected Hamiltonian in the first order approximation is 
given by 

H = E M | $n [n] )($ || w | + J x P {l J2Si,iSi,i+iP\\ + 

[n] i,l 

JdP\\ ^(SuS i+ ij + i + S i+ ijSid + i)P\f21) 

i,l 

which may simply be written as 

[n] il 

+ Sj+ijSy+i), (22) 

il 

where S™ ! ^ m! = (0 n , \S it i \4> mi )- 

The matrix elements for the second terms may be writ- 
ten 



(23) 



I 

ment to be non zero is that (4>nn \4>mi 2 ) — ^nn,mi2 except 
when Z12 = I or I + 1. Thus, 



($ U[n] | J ffj.Qyjr x |$|| [n , ] )= J2 

, [m] 



|Si, Z S M+ l ($|| [ m ] |S M S M+ 1 1$|| 



E 



II [0] 



En 
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E [rn] - E[0] - ( e m - EOi) + ■•• + fan, ~ e O, ) + (e TO; + 1 - £0 ; + 1 ) + ••• + {(-n L - (Q L )- (24) 

I 



FIG. 2: sketch of the superblock in the step 1 

In the Eq. 1241 above, the dominant terms will come 
from the differences involving the indices mi and m;+i 
because the others terms come from the state used to 
generate P|| and are thus of lower energies. Up to the 
second order, the effective one-dimensional Hamiltonian, 
which is written here without the frustration term, is 

H ~ Z! %«] 1*11 W>(*ll W I + J ^ Z - 

[n] i 

4^Sp)sS + ... (25) 
I 

where the chain-spin operators on the chain Z are Si = 
(Sh,S 2 ,,...S w ) and S<j° = (Sg>, §£>,...§$), L is the 
chain length. The matrix elements of the second order 
local spin operators are 

Sf nin > = Y ^ ^ — . (26) 

(2) 

One can note that this expression of S a is not exact, it 
has been simplified to avoid long-range coupling between 
the chains. The effective ID Hamiltonian H is also stud- 
ied using DMRG. 

IV. ALGORITHMIC DETAILS 

The algorithm of the method will now be described 
below. It consists of two DMRG steps separated by an 
intermediate stage in which a simple block decimation is 
made. 



A. Step 1 

The first step of the method is the usual DMRG 
method for a single chain. The chain is divided into four 
blocks, and the two internal blocks are made of a single 
site each. In the calculations, ms\ = 16 — 192 states are 
kept in the two external blocks. In most cases, the ini- 
tial iteration starts with a chain having the largest size 
before truncation, for instance L = 16 when ms\ = 128 
states are kept. This way, a high accuracy is obtained 
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FIG. 3: sketch of the transformation of the two external 
blocks of length L/2 into a single block of length L which 
is used as the building unit in step 2 

even when the infinite system method is used. During 
this step, the local spin operators St on each site i of the 
chain are stored and longitudinal spin-spin correlations 
C(i,r) = (SjSi_|_ r ) are also computed and stored. As 
discussed by Caron and Bourbonnais 2 ^, open boundary 
conditions (OBC) which are used here introduce spurious 
behavior at the edges of the chain. It is therefore better 
to chose the origin i in C(i,r) in the middle of chain. It 
is crucial during this step to target more just than the 
S z = sector in order to obtain a correct low-energy 
Hamiltonian. In addition to S z = 0, S z = ±1, ±2 were 
targeted in this study. 



B. Block transformation 

An intermediate stage of the algorithm is a decima- 
tion process as in the old block renormalization group 
method^iiS. In this process, the two external blocks hav- 
ing L/2 sites each are reduced to a single block with L 
sites. During this step, the ms\ x ms\ states describing 
the chain are reduced to ms 2 lowest states of the chain. 
As noted in Refi, since the block transformation is used 
only one time, the problem of the propagation of spurious 
boundary effects^ , is not present. All the local spin op- 
erators and spin-spin correlation functions are expressed 
in the basis of the ms 2 states. 



C. Step 2 

The second step consists of applying the ID DMRG 
method using the chains obtained at the end of the pre- 
vious step as the building blocks. This step is indeed 
identical to the first step, except for the dimension of the 
local spin operators. The central block is the chain from 
the previous step and thus has the dimension ms 2 x ms-i. 
Typically, msi = 16 — 96 and for the two external blocks, 
roughly the same number of states is kept. If four blocks 
were taken as in the first step, the dimension of the su- 
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FIG. 4: Sketch of the superblock in step 2 

perblock would be ms\ which can become rapidly im- 
pratical. To ease the computations, three blocks instead 
of four are mostly used during this step. As will be seen 
below, it is important during this step to check that, for 
a given value of the couplings J± and J^, enough states 
are kept such that a valid computation is made, i.e., 
that the truncated Hamiltonian generated for the single 
chain is accurate enough, for the ground state and for 
the low lying states, to be used as a building block for 
the 2D lattice. One can easily see that for a fixed L and 
J±.d <C A a (L), A a (L) is the finite size spin gap, and the 
interchain matrix elements will be negligible. The sys- 
tem will behave as a collection of free chains even if Jx,d 
is turned on. Now if Jj_,d ~ 8E(L), where 6E(L) is the 
width of the retained states, the matrix elements of the 
states having higher energy, which have been truncated 
out, have a non-negligible contribution. 

D. Algorithm 

The algorithm is summarized below. 

• 1. Build the low energy Hamiltonian for a single 
chain B ■ -B by using the ID DMRG algorithm of 
Refi. Store the spin operator S, on each site and 
the correlation function C(i,r). 

• 2. When the block B size is L/2, apply the block 
method to merge the two external blocks into a sin- 
gle block defined by the msi states kept. Express 
all the spin operators and correlation functions in 
the basis of the ms2 states. Check if for the num- 
ber of states kept, the transverse couplings satisfy 
A a (L) < Jj_^ <C 5E(L). If this condition is not 
satisfied, increase ms2- 

• 3. Start a second ID DMRG simulation identical 
to the first one except that the central block is now 



J 


TllS 2 — 16 


ms2 = 24 


771 S 2 — 32 


QMC 


0.00 


-0.42848 


-0.42851 


-0.42851 


-0.42849(2) 


0.05 


-0.42900 


-0.42907 


-0.42909 


-0.42926(2) 


0.10 


-0.43058 


-0.43078 


-0.43090 


-0.43147(2) 


0.15 


-0.43312 


-0.43361 


-0.43387 


-0.43530(2) 


0.20 


-0.43642 


-0.44733 


-0.43780 


-0.44064(2) 


0.25 


-0.44028 


-0.44174 


-0.44247 


-0.44727(2) 



TABLE I: DMRG ground state energies for 12 x 12 lattices 
for msi = m,S2 = 16, 24, and 32 versus QMC. 



a single chain instead of a site, and the exchange 
coupling is Jx,d instead of 1.0. 

V. RESULTS WITH FOUR BLOCKS IN STEP 2 

In this part, the DMRG results are compared to the 
stochastic series expansion (SSE) QMC results. The 
SSE-QMC method^ is so far the most reliable tech- 
nique for the study of quantum spin systems. It has been 
used to study weakly coupled quantum spin chains^. It 
will thus be very instructive to see how well the DMRG 
method compares to the SSE-QMC. 



A. First order ground-state energies 

In Table (HJ, the ground state energy per site for 12 x 
12 systems for m — 16, 24 and 32 is shown. In this 
calculation, four blocks were used in the second DMRG 
step. The agreement with the SSE-QMC results is good 
for small g. The DMRG energies are higher than those of 
the QMC for all transverse couplings studied. As ms2 is 
increased, the difference between the DMRG and QMC 
energies decreases. This was expected since the current 
method as the original DMRG procedure is variational. 

The band-width of the states kept is AE — 1.132, 
1.290, 1.466 when ms2 — 16, 24 and 32 respectively. 
The target states during the first DMRG step were the 
lowest states of the spin sectors with S z = 0, ±1 ±2. 
The lowest states of higher spin sectors have energies 
which are higher than the highest state kept in lower 
spin sectors, therefore they were not targeted. The fact 
that the DMRG results compare well with the QMC ones 
even at intermediate couplings reveals that for the spin 
chain, reliable calculations can be made for values of 
SE(L)/J± w 5. But as expected for higher values of 
J±, the condition 5E(L)/J± ^ 1 is no longer fulfilled. 
This means the Hilbert space is too severely truncated. 



B. Second order ground-state energies 

Table [H] displays second order ground-state energies 
for a 12 x 12 system which are compared with QMC. The 
agreement is systematically better than in the first order 
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" _L 


777, So — 16 


ms2 = 24 


tiis 2 — 32 


QMC 


0.00 


-0.42848 


-0.42851 


-0.42851 


-0.42849(2) 


0.05 


-0.42901 


-0.42909 


-0.42910 


-0.42926(2) 


0.10 


-0.43063 


-0.43083 


-0.43094 


-0.43147(2) 


0.15 


-0.43322 


-0.43369 


-0.43396 


-0.43530(2) 


0.20 


-0.43661 


-0.43746 


-0.43797 


-0.44064(2) 


0.25 


-0.44055 


-0.44192 


-0.44271 


-0.44727(2) 



TABLE II: DMRG ground state energies for 12 x 12 lat- 
tices for second order DMRG compared with QMC for msi = 
ms 2 = 16, 24 and 32. 



case for all values of J± studied. But the improvment 
is small. This is because as discussed above, the DMRG 
does not provide the full ID spectrum. Only the states 
kept to form the reduced superblock are used in the per- 
turbative expansion. When ms2 = 32, this is merely 924 
states i.e., a tiny fraction of of the 2 144 states which form 
the full Hilbert space of the 12 x 12 lattice. Another 
reason for this modest improvment is the fact that the 
DMRG energies are variational. The high lying energies 
which are used to generate second order terms are ob- 
tained with less accuracy than the states kept in the first 
order. Indeed, this does not mean that the matrix ex- 
pansion presented above is not efficient. It has been used 
in the simple case of the Mathieu equation for which the 
full spectrum of unperturbed Hamiltonian is available^. 
The convergence of the matrix method is quite impres- 
sive. Thus it seems that the best way to use the matrix 
Kato-Block expansion when the DMRG is used to obtain 
the unperturbed spectrum is to restrict oneself to the first 
order and keep msi as large as possible. However, larger 
values of msi are unpratical when four blocks are used to 
form the superblock in the second step. For this reason 
from now, only three blocks will be used to generate the 
superblock in the second step. 



VI. FIRST ORDER RESULTS WITH THREE 
BLOCKS IN STEP 2 

When three blocks are used, the superblock size in the 
second step is divided by ms2 relative to the case of four 
blocks. This significantly reduces the amount of required 
CPU for a given value of ms2- But this is not with- 
out problems. It was noted thal£, when three blocks are 
used, the coupling between blocks may incorrectly sets 
in leading to a poor performance of the method even 
if the truncation errors are small. The remedy against 
this problem is to target more than one state so that the 
interblock mixture is performed correctly. However, tar- 
geting many superblock states lower the accuracy on the 
ground state. For this reason, only the ground state was 
targeted. The truncation errors were in general smaller 
than 10~ 8 for ms2 varying from 16 to 96 for different 
lattice size. But as said above, this does not give any 
indication about the accuracy of the second step of the 



J 


777, S 2 — 16 


777, 5 2 — 32 


77i So — 64 


QMC 


0.00 


-0.42187 


-0.42187 


-0.42187 


-0.42186(2) 


0.05 


-0.42239 


-0.42244 


-0.42247 


-0.42246(2) 


0.10 


-0.42402 


-0.42421 


-0.42439 


-0.42444(2) 


0.15 


-0.42670 


-0.42722 


-0.42762 


-0.42771(2) 


0.20 


-0.43032 


-0.43144 


-0.43219 


-0.43239(2) 


0.25 


-0.43470 


-0.43673 


-0.43799 


-0.43843(2) 


TABLE III: DMRG 


ground state 


energies for 8x9 lattices 


for msi 


= 16, and m,S2 = 16, 32, and 64 versus 


QMC. 


Jl 


ms2 = 32 


ms2 = 64 


ms 2 = 80 


QMC 


0.00 


-0.42851 


-0.42851 


-0.42851 


-0.42850(1) 


0.05 


-0.42910 


-0.42918 


-0.42919 


-0.42922(1) 


0.10 


-0.43094 


-0.43124 


-0.43131 


-0.43150(1) 


0.15 


-0.43396 


-0.43468 


-0.43483 


-0.43537(1) 


0.20 


-0.43796 


-0.43928 


-0.43956 


-0.44075(1) 


0.25 


-0.44268 


-0.44476 


-0.44521 


-0.44744(1) 



TABLE IV: DMRG ground state energies for 12 x 13 lattices 
for msi — ms 2 = 32, ms2 = 64 and ms2 = 80 (msi = 64 for 
both) versus QMC. 



method. The QMC results are thus taken as the reference 
to gauge the DMRG results. 



A. Ground-state energies 

For small sizes and weak couplings, the differences be- 
tween the DMRG and QMC ground state energies are 
very small. For instance for the 8x9 lattice shown in 
Table( llllj) . for J± = 0.05, the difference between the 
two methods is only 0.00016 for ms2 = 16. The two re- 
sults are within QMC error when ms2 is increased to 64. 
As expected, increasing the coupling tends to decrease 
the accuracy because the ratio 6E(L)/J± is reduced. In- 
creasing the lattice size has the same effect on this ratio 
because 8E(L) is smaller for larger lattices for a fixed 
ms2 (Tabic ( ITTT1 IIVI Ivl) ) . One may note that by keeping 
a larger number of states than in the case of four blocks, 
the accuracy has increased in all cases. 



Lx L+l 


DMRG 


QMC 


8x9 


-0.42440 


-0.42444(2) 


12 x 13 


-0.43124 


-0.43150(2) 


16 x 17 


-0.43481 


-0.43529(1) 



TABLE V: DMRG ground state energies for various lattices 
for ms2 — 80 and J± = 0.1. 
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I 


DMRG (1) 


QMC (1) 


DMRG (t) 


QMC (t) 


1 


-0.14595 


-0.14931(1) 


-0.02047 


-0.02209(1) 


2 


0.06072 


0.05904(1) 


0.00561 


0.00525(1) 


3 


-0.04799 


-0.05173(1) 


-0.00191 


-0.00164 


4 


0.03340 


0.03537(1) 


0.00066 


0.00055 


5 






-0.00023 


-0.00019 



TABLE VI: DMRG versus QMC longitudinal (1) Cy (7, 7, r) 
and transverse C±(7,7,r) spin-spin correlations for a 12 x 13 
lattice for J± — 0.1, Jd = 0. 



B. Ground-state correlation functions 

It is not possible to keep track of all spin-spin corre- 
lations when large systems are studied because of CPU 
and memory limitations. The behavior of spin-spin cor- 
relations is thus studied along one chain in the direction 
parallel to the chains and one chain in the direction per- 
pendicular to the chains. These correlation functions are 
respectively given below: 

q|(» > / > r) = i<S i ,,.S i+r ,j). (27) 
C ± (*,l,r) = i<S i) ,.S<, I+r ) (28) 

It is particularly difficult to obtain the large r behav- 
ior of the correlation functions because of a number of 
factors that complicate such an analysis. At the level of 
a single chain, the long distance behavior of C(i, r) is al- 
ready complicated by logarithmic corrections. Although 
highly accurate data can be obtained in ID from QMC 34 
or DMRG™, the two studies disagree on the exact form 
of the logarithmic corrections. Furthermore when open 
boundary conditions (OBC) are used instead of periodic, 
the spin-spin correlation functions show strong odd-even 
alternations 2,29 . This is because the ground state may be 
regarded as a resonant state between a state with strong 
bonds on even links and weak bonds on odd links, and 
a state with weak bonds on even links and strong bonds 
on odd links. Another difficulty with OBC is that the 
translational invariancc of the chain is broken, and the 
value of Cu(i,l,r) depends on the position of the site 
chosen as the origin on the lattice. It was shownSS that 
the closer the origin is to the edge of the lattice, the 
higher are the spurious effects introduced by the OBC. 
All these facts render the direct detection of long range 
order in the transverse direction, for which the spin-spin 
correlations are very small, impossible to achieve with the 
present calculation for which the magnitude of C±(i, I, r) 
for large r is close to the accuracy on the eigenvalues dur- 
ing each iteration. An alternative way is to look at the 
Cj|(i,Z,r), because the existence of long range order in 
the longitudinal direction is an indication that the order 
is two-dimensional. 

In order to observe the correct long-range behavior, 
one must first reduce the influence of the spurious effects 



I 


DMRG (1) 


QMC (1) 


DMRG (t) 


QMC (t) 


1 


-0.14640 


-0.14619(1) 


-0.02116 


-0.02533(1) 


2 


0.06059 


0.06130(1) 


0.00726 


0.00854(1) 


3 


-0 04-875 


-0 odQssn i 


-0 00390 


-0 00399 


4 


0.03422 


0.03537(1) 


0.00147 


0.00201 


5 


-0.02866 


-0.02990(1) 


-0.00078 


-0.00105 


6 


0.02251 


0.02363(1) 


0.00030 


0.00056 


7 






-0.00013 


-0.00030 


8 






0.00006 


0.00015 



TABLE VII: DMRG versus QMC longitudinal (1) C\\ (9, 9, r) 
and transverse Cj_(9, 9,r) spin-spin correlations for a 16 x 17 
lattice for J± — 0.1, Jd = 0. 

generated by the application of the OBC. Furthemore 
to simplify the analysis, the eventual logarithmic correc- 
tions will not be considered here. In order to avoid the 
odd-even alternation, Cm (i, I, r) and Cj_(i, I, r) were aver- 
aged in the period of these alternations. This was done 
by computing (S^jSi+^j) at two different origins. The 
spin Su is taken as the origin of a strong link or as the 
origin of a weak link. The actual correlation function is 
then 

C|| (t, I, r) = 0.5(C|| (i, I, r) + C\\ {% + 1, I, r)). (29) 
And for C±(i, l,r), 

C±(i ) l,r) = 0.5(C±(i,l,r) + C±(i ) l-l,r)). (30) 

The averaged correlations Cm (i, I, r) and C±(i, I, r) are 
shown in Tables IVII IVIII and IVIIII for, respectively, 
12 x 13, 16 x 17 and 32 x 33 lattices. The origins (£, I) 
of the correlation functions were chosen at the middle 
of the chain in order to minimize the end effects, (i, I) 
was equal to (7, 7), (9, 9) and (17, 17) respectively for the 
12 x 13, 16 x 17 and 32 x 33 lattices. For the 16 x 17 lat- 
tice, ms\ = 128 states were kept during the first DMRG 
step and ms2 = 64 states were kept during the second 
DMRG step. The comparison with QMC is quite good 
in the longitudinal direction but less good in the trans- 
verse direction when the lattice size gets large. For the 
32 x 33 lattice, ms\ and ms2 were respectively increased 
to 160 and 80. As for the case of the 16 x 17 lattice the 
agreement was quite good for Cm (i, I, r) and less good for 
C±(i, l,r). The reasons for the differences are not easy to 
analyze. Although very small truncation errors p m (for 
instance, p m < 1 x 10~ 7 for msi = 128 and ms2 = 64)) 
are obtained in the DMRG, there is no obvious relation 
between these truncation errors and the errors on the 
measurements. Furthermore, the effects of higher order 
terms in the perturbation series have not be analyzed for 
the case of three blocks. Since more states are kept when 
three blocks are used, the contribution of second order 
terms is likely larger than the one found above for four 
blocks. 



1 


DMRG (1) 


QMC (1) 


DMRG (t) 


QMC (t) 


1 


-0.14694 


-0.14636(3) 


-0.01846 


-0.02952(2) 


2 


0.06042 


0.06151(3) 


0.00969 


0.01465(2) 


3 


-0.04908 


-0.05066(2) 


-0.00623 


-0.01057(3) 


4 


0.03402 


0.03640(2) 


0.00416 


0.00821(1) 


5 


-0.02949 


-0.03229(3) 


-0.00281 


-0.00662(1) 


6 


0.02366 


0.02682(3) 


0.00190 


0.00545(1) 


7 


-0.02108 


-0.02450(2) 


-0.00128 


-0.00453(2) 


8 


0.01820 


0.02163(3) 


0.00086 


0.00379(2) 


9 


-0.01643 


-0.01990(2) 


-0.00059 


-0.00321(2) 


10 


0.01474 


0.01806(2) 


0.00040 


0.00270(2) 


11 


-0.01327 


-0.01646(2) 


-0.00028 


-0.00228(2) 


12 


01 202 


01 !iOS,(2) 


0001 8 


001 Q4f2l 


13 


-0.01045 


-0.01309(1) 


-0.00012 


-0.00164(2) 


14 


0.00914 


0.01164(2) 


0.00008 


0.00137(2) 


15 






-0.00005 


-0.00115(1) 


16 






0.00003 


0.00088(1) 



TABLE VIII: DMRG versus QMC longitudinal (1) 
Cm (17, 17, r), and transverse (t) Cj_(17, 17, r) spin-spin 
correlations for a 32 x 33 lattice for J±_ = 0.1, Jd = 0. 



VII. GROUND-STATE PROPERTIES IN 
PRESENCE OF FRUSTRATION 

The DMRG method has shown an overall good agree- 
ment with QMC for weak couplings and not too large 
sizes. The method is well controlled and can system- 
atically be improved by increasing ms\ and msi. The 
advantage of the DMRG over QMC is that it is very 
flexible and can be applied to frustrated systems. A sit- 
uation where the QMC is known to fail. In this section 
a diagonal Jd exchange coupling is included. It has the 
effect of introducing a competition between interchain 
AFM correlations along the rows and AFM correlations 
along the diagonals. 



A. Ground-state energies 

Although the result on the ground state energy can not 
provide information about a possible long-range order, it 
can be helpful to see if the perturbation is relevant or 
not. Fig. [3 shows the binding energy per chain Eb = 
E (L)-E (Lx(L+l))/(L+l), where E (L) is the ground 
state energy for a single chain and Eq(Lx (L + 1)) is the 
ground state energy for an L x (L + 1) lattice. J± is set to 
0.1. Eb first decreases as Jd/J± is increased. It reaches 
a minimum at Jd ~ 0.5Jj_. At the minimum point, the 
binding energy nearly vanishes, Eb ~ 0.0015 which is 
roughly two orders of magnitude smaller than its value 
for Jd — 0. As Jd/J± is further increased, Eb starts 
increasing sharply. This behavior suggests the existence 
of three regimes for for the action of small perturbations 
on the single chain, two stable phases separated by a 
transition region. The first regime, which occurs when 
Jd 0.5 J±, is a Neel state as is already known from 
QMC studies 14 . This will be confirmed below by the 




P 

FIG. 5: The binding energy Eb with respect to single chain 
for a 32 x 33 lattice as a function of p = Jd/J±, J± = 0.1. 
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FIG. 6: The binding energy Eb with respect to single chain 
for a 32 x I lattice as a function of I for Jd/J± = (circles), 
Jd/J± = 0.5 (squares), Jd/J± = 1-2 (diamonds), J± = 0.1. 



analysis of the DMRG correlation functions. The second 
regime is when Jd ~ 0.5 Jx- The perturbation seems 
to be irrelevant, J± and Jd cancel each other so that 
there is almost no gain in energy by applying the two 
perturbations simultaneously. In the third regime, when 
Jd ^ 0.5Jj_, the ground state is also magnetic with a 
collinear order, an alternate arrangement of transverse 
up and down ferromagnetic chains (see Fig.^). 

The above analysis is further supported by observing 
the evolution of the binding energy Eb(L x I) as a func- 
tion of the number of chains in the lattice (Fig. EJ. It 
clearly shows that when Jd ~ 0.5Jj_, the binding energy 
is nearly independent of the number of chains and re- 
mains very close to that of the single chain. Hence it 
seems that at the point Jd ~ 0.5Jj_, the ground state is 
made of independent chains as for J± = Jd = 0. This be- 
havior is analogous to the domino model studied by Vil- 
lain and coworkers^, where a disordered ground state, 
made of independent chains for a particular value of the 
transverse coupling, was found. 



0.4 r 1 1 1 1 1 1 1 1 1 1 r 




FIG. 7: The ground state transverse structure factor for a 
32 x 33 lattice for Jd/J± — (filled circles), 0.2 (open circles), 
0.4 (filled squares), 0.55 (open squares), 0.8 (filled diamonds), 
1.0 (open diamonds), 1.2 (open triangles). 



B. Ground-state correlation functions 

The behavior of the correlation functions is consistent 
with the existence of the three regimes found for the 
ground state energy. As expected, spin-spin correlations 
along the chains remain antiferromagnetic. The change 
of regimes will be detected by analyzing spin-spin corre- 
lations along the transverse direction. Fig. [7| shows the 
transverse magnetic structure factor S±(k±), 



L/2 

S±(k±_)= C ± (17, 17, r) cos k±r (31) 

kj_=l 

where k± is a wave number in the transverse direction. It 
also shows the three regimes discussed above. When Jd < 
Q.5Jx, Sj_{k±) has a maximum at fcj_ — n. The spin- 
spin correlations along the tranverse direction are AFM 
as for the longitudinal direction. For Jd « 0.5Jj_, S± is 
structureless, a fact which is consistent with disconnected 
chains. When Jd 0.5 Jj_, S± has a maximum at k±_ = 0, 
and the correlations in the transverse direction are now 
ferromagnetic. This is the collinear magnetic state shown 
in Fig.[T] 

The bond-strength Cj_(17, 17, 1), computed in a 32x33 
lattice is shown in Fig. |SJ It also shows that the chains 
seem to be disconnected when Jd ~ 0.5Jj_. Starting 
from J d = 0, for which Cj. (17, 17,1) = -0.01846, its 
absolute value first slowly decreases. Then, when Jd is 
in the vicinity of Jd c , the absolute value of Cj_(17, 17, 1) 
sharply decreases and become very small; Cx(17, 17, 1) = 
—0.000799 when Jd = 0.5Jj_. As soon as Jd exceeds 
Jd C : Cj. (17, 17,1) becomes ferromagnetic and starts to 
increase sharply. It later saturates when one is far enough 
from the critical point. 
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FIG. 8: The bond-strength C(l) = C ± (17, 17, 1) as a function 
of p = J d /J±, J± = 0.1. 



VIII. LONG-RANGE ORDER IN THE GROUND 
STATE 

The analysis made in the preceeding section indicates 
regions of dominant Neel or collinear spin-spin correla- 
tions or of a possibly disordered ground state at the tran- 
sition point. But it does not tell if long range order is 
truly established. For this, it is necessary to look at the 
long-range behavior of the correlation functions. 

The spurious effects due to the breaking of the trans- 
lational symmetry, a consequence of the OBC, may be 
reduced by using a filter which smooths the action of the 
sites near the edges. In the results shown below in Fig 
El El El this was done as follows: Cn(i, Z,r) was first 
examined for a single chain for which the long distance 
behavior is known. Roughly, Cn(i,l,r) cx 1/r if logarith- 
mic corrections are neglected. It was found that if the 
origin is taken at the middle of the chain, the 1/r behav- 
ior is roughly satisfied for d < r < L/2 — d with d w 8. 
The second inequality is due to edge effects. As a conse- 
quence, relatively large values of L are necessary in order 
to observe the long range behavior, and lattices of up to 
64 x 65 were studied. The problem with such large lattices 
is that the energy width SE(L) shrinks with increasing 
L and the condition A a (L) ~ J±,d -C 5E(L) may not be 
fulfilled. For L — 64, ms\ — 192 and ms-i = 80 states 
were kept. For these values, 8E(L — 64) w 0.5, which 
means SE(L = 6i)/J± : d ~ 5 provided that J± : d ^5 0.16. 

The first question which needs to be addressed is to 
know whether the DMRG can detect an eventual long 
range order. Comparisons with QMC for L — 32 show 
that, the DMRG correlation in the transverse direc- 
tion decays faster. This effect is expected to be larger 
on longer chains. But, despite this shortcoming of the 
method, one can still detect possible occurrence of long 
range order. If one considers the central chain in the 2D 
lattice, Ch(L/2 + 1, L/2 + 1, r) is modified from that of an 
isolated chain because of the effective magnetic field cre- 
ated on it by the rest of lattice. Although this effective 
field is somewhat under evaluated by the DMRG because 
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FIG. 9: The ground state correlation function Cu (33, 33, r) 
for the 64 x 65 lattice for J± = (circles) and Jj_ = 0.16 
(squares). Jd = in both cases. The filled and open symbols 
correspond to odd and even distances respectively. 



the transverse correlations are underevaluated, it can still 
be strong enough to lead to an ordered phase. This inter- 
pretation is related to the chain mean- field approach 36 ; 
the essential point is that, here, no assumption about 
long-range order is made a priori. From this, one can 
see that if the DMRG method leads to a finite order 
parameter, it is necessarily genuine. Fig. compares 
for the correlation function Cm (33, 33, r) for J± = 0, 
Jd = with J± = 0.16, Jd = 0. In the first case when 
both transverse couplings are absent, Cm (33, 33, r) oc 1/r. 
The DMRG data still show an odd-even alternation, so 
fits must be perfomed for odd and even distances sep- 
arately. The best least square fits to the data gave 
Cm (33, 33, r — > oo) w 0.0001. This is consistent with an 
absence of a long-range order for an isolated chain. But 
in the case J± = 0.16 and Jd = 0, a fit to the data shows 
that C||(33,33,r) tends to C||(33,33,r -► oo) « 0.0036 
The existence of long-range Neel order for J± = 0.16 is 
consistent with previous studie d 12 ' 14 -. Adding J± alone 
seems to lead to long-range order. This has been re- 
cently shown in Refi 1 ^ where values of J± down to 0.02 
were investigated. It is of course impossible to show from 
a numerical investigation whether any small value of Jj_ 
will lead to an ordered state or there may be a disordered 
state for very small values of J±. In view of current nu- 
merical results, the former hypothesis is more convincing. 

The above discussion suggests that a frustration Jd 
must be added in order to thwart the Neel state which 
results from the action of J±. J± will now be set to 
0.16 and Jd varied. For Jd — 0.08, the value at which 
the analysis of smaller chains suggested that the ground 
state is made of disconnected chains, C\\ (33, 33, r) is com- 
pared to the same quantity for a single chain in Fig. 1101 
Fits to the data show that the behavior of C|| (33, 33, r) 
is quite similar to that of a single chain. Clearly, for 
these values of the transverse couplings there is no long- 
range order in the ground state, and Cm (33, 33, r) seems 
to indicate that the ground state is made of a set of in- 
dependent chains. It is important to emphasize that this 
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FIG. 10: The ground state correlation function Cu (33, 33, r) 
for the 64 x 65 lattice for J± — 0., Jd = 0. (circles) and 
J± — 0.16, Jd = 0.08 (squares). The filled and open symbols 
con 
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FIG. 11: The ground state correlation function Cm (33, 33, r) 
for the 64 x 65 lattice for J± — 0., Jd = 0. (circles) and 
J± — 0.16, Jd = 0.16 (squares). The filled and open symbols 
correspond to odd and even distances respectively. 

result does not mean that the finite temperature effects 
are also trivial. The present situation could be similar to 
the so-called domino model first introduced by Andrei 
and later studied by Villain and coworkers^ or to the 
crossed-chains quantum spin models^. In the domino 
model, it was found that the ground state was made of 
disconnected chains but there was a long-range order at 
finite temperature. Indeed, the Mermin- Wagner theorem 
prohibits long-range order at finite temperatures for the 
2D Heisenberg model. The finite temperature behavior 
in this case will thus be different. 

The disconnected chain ground state is in contradic- 
tion with a recent study by Nersesyan and Tsvelik 15 . 
These authors argued, using bosonization, that when 
Jd/ J±. — 0.5, only the staggered part of the interchain 
part of the Hamiltonian vanishes. There remains a uni- 
form part which is relevant and leads to two-dimensional 
spin liquid with a spin gap, A CT oc exp(— fjj-), where 
v a is the spin velocity. The low energy excitations are 
argued to be unconfined spinons. The apparent contra- 
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diction between this conclusion and the numerical data 
above could be that the binding energies of the 2D spin 
liquid are very small, indeed J± = 0.1 corresponds to 
A CT « 1.0 x 10 -11 . Such a small energy can obviously 
not be detected by a numerical method. A way to avoid 
this small energy scale is to raise J±, This possibility is 
currently being investigated. 

Finally, the collinear magnetic long range order is also 
confirmed by the analysis of (33, 33, r-) . In Fig. ^2 h 
is shown that for Jd — 0.1, C\\ (33, 33, r) converges even 
faster than for the Neel state above. The value of the 
extrapolated correlation is C*||(33, 33, r — > oo) w 0.0056. 

IX. CONCLUSIONS 

In this paper, a new renormalization group method 
for weakly coupled chains was presented. It is based on 
solving numerically the model Hamiltonian in two ID 
steps using the DMRG. During the first step, a low en- 
ergy Hamiltonian for a single chain is obtained using the 
ID DMRG. The original problem is then formulated as 
a perturbative expansion around the DMRG low energy 
Hamiltonian obtained during the first step. This per- 
turbative expansion is a ID problem which can also be 
solved by the DMRG. 

The first and second order approximations were stud- 
ied for weakly coupled Heisenberg chains with and with- 
out frustration. The results were compared to the QMC 
and showed good agreement for small systems and small 
transverse couplings. It was shown that, starting from 



the disordered ID chain, the method can predict long- 
range order when it exists, a test generally failed by con- 
ventional perturbative methods. Calculations performed 
in the presence of frustration indicate an absence of a 
genuinely 2D spin liquid state. Instead, the frustration 
drives the Nccl ground state to a collinear magnetic state. 
At the transition point, both ground-state energy and 
spin-spin correlation functions show a disordered ground 
state. The precise nature of this disordered ground state 
is currently under investigation. 

The above results are very encouraging and indicate 
that the DMRG may become a very useful tool for the 
study of highly anisotropic 2D systems in the future. 
The method is only in its early stages, and some im- 
portant improvements of the method are currently un- 
derway. These are the investigation of the role of cluster 
corrections, i.e., the starting point in the first step will 
be two-leg or three-leg ladders instead of a single chain; 
the use of exact diagonalization during the first step in- 
stead of DMRG. These improvements are likely to lead to 
better results for spin-spin correlations in the transverse 
direction. Extensions of the method to thermodynamic 
spin systems or fermionic models will also be made in the 
near future. 
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